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Abstract 

We describe a new method for finding analytic solutions to some initial-boundary problems for partial differ- 
ential equations with constant coefficients. The method is based on expanding the denominator of the Laplace 
transformed Green's function of the problem into a convergent geometric series. If the denominator is a linear 
combination of exponents with real powers one obtains a closed form solution as a sum with finite but time 
dependent number of terms. We call it a D'Alembert sum. This representation is computationally most effec- 
tive for small evolution times, but it remains valid even when the system of eigenmodes is incomplete and the 
eigenmode expansion is unavailable. Moreover, it simplifies in such cases. 

In vibratory problems D'Alembert sums represent superpositions of original and partially reflected traveling 
waves. They generalize the D'Alembert type formulas for the wave equation, and reduce to them when original 
waves can undergo only finitely many reflections in the entire course of evolution. The method is applied to 
vibrations of a bar with dampers at each end and at some internal point. The results are illustrated by computer 
simulations and comparisons to modal and FEM solutions. 

Keywords: wave equation, viscous boundary conditions, superposition of traveling waves, D'Alembert's 
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1 Introduction 



We describe a new method for finding analytic solutions to initial-boundary problems for partial differential equa- 
tions with constant coefficients, and apply it to vibrations of a bar with two dampers at the ends and one at an 
internal point. We start, as in other approaches, by taking the Laplace transform with respect to time and find 
the Green's function of the resulting boundary eigenvalue problem. The general idea is to expand it into a series 
over functions with simpler dependence on the Laplace parameter, and then to invert this series termwise. In the 
modal approach, which is most commonly used pQ, these simpler functions are the partial fractions. Termwise 
inversion leads to an expansion over the eigenmodes (standing waves) of the problem, which often have explicit 
Laplace inverses. Our expansion also uses functions with explicit inverses, but its terms are weighted by exponents 
with negative real powers. The latter invert into time shifted Heaviside functions and insure that only finitely many 
terms of the expansion contribute to the solution at any given time. As a result, for any finite time our method 
produces closed form solutions, finite sums that we call D'Alembert sums. Analytically, the matter reduces to ex- 
panding the denominator of the boundary eigenvalue problem's Green's function into a convergent geometric series. 
One essential restriction is that this denominator must be a linear combination of exponents with real powers. 

While this approach to Laplace inversion is common for ordinary differential equations, e.g. in signal processing 
2 j Sec 3.11], it docs not seem to be popular for partial differential equations. One reason is that the modal expansion 
is more universal. Another is that one can achieve a resonable approximation for arbitrary times with a limited 
number of eigenmodes, while the number of terms in a D'Alembert sum may grow fast with time. However, when 
the method works it produces an exact solution, and more importantly, it works particularly well when the modal 
approach breaks down. Recall that a system of eigenmodes may be incomplete, i.e. not all functions can be 
approximated by linear combinations of eigenmodes, for some values of problem parameters called critical [3 . In 
the critical cases modal expansion does not exist, but a D'Alembert sum not only exists but takes a particularly 
simple form. 

Rather than to begin with an abstract exposition, we first show how the method works in simple examples 
involving a vibrating bar with dampers attached at the boundaries and at some internal point. Particular instances 
of this initial-boundary problem served as the main motivation for this work. Aside from being a perfect illustration, 
this problem provides some helpful physical insight into how and why the method works. To be more specific, let 
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us describe the problem in more detail. 



u(x, t) 




Figure 1: A bar with viscous ends and internal damper. 



Figure [T] depicts a bar of length L, free to move horizontally, suspended by two dampers at each end and by one 
at the distance a from its left end. Symbols p 1 A and E represent the density of the bar, the constant cross-sectional 
area and its modulus of elasticity respectively, the wave speed along the bar is denoted c := (E/p) 1 / 2 . Let c\, C2 and 
C3 be the damping coefficients of the left, right and internal dampers respectively, we set h\ :— j^ci, /12 := -§^C2 
and /13 := 2§a c 3 ^ ne ex t ra 1/2 simplifies some formulas). These hi along with a/L are the dimensionless parameters 
that determine qualitative behavior of the system. Since in our case the bar can move rigidly, just as in the problem 
with free ends, we write the equations of motion in the absolute frame that remains at rest at all times. At t = 
the left end of the bar is assumed to be at the origin, and u(x, t) denotes the displacement of the point with initial 
coordinate x at time t in the absolute frame, see FigfTJ 

The system is described by a modified wave equation 

u t t(x,t) + 2h 3 c5(x - a)u t {x,i) - c 2 u xx (x,t) =p(x,t), (1) 

with the boundary conditions 

« B (0,t)-^u t (0 I t) = and u x (L, t) + -± u t {L, t) = 0. (2) 

Here p{x, t) is the external force per unit mass, and the subscripts x , t denote partial derivatives with respect to space 
and time. Given u(x, t), the solution in the frame that moves along with the left end of the bar is u(x,t) — u(0,t). 
From the modal point of view the problem was analyzed in [4]. As long as hi,fi2 ^ ±1 it has a countable discrete 
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spectrum with a complete set of eigenmodes. When one or both hi = — 1 the problem is superunstable and solution 
exists only for special initial data and for finite time only, cf. [S]. In the remaining critical cases, hi = 1 and/or 
I12 = 1, one or both boundaries become transparent, i.e. no energy is reflected back by them. As a result, the 
standing waves can not form along the entire length of the bar, explaining why the system of eigenmodes becomes 
incomplete and modal expansion becomes impossible. When the standing waves can not form even on parts of the 
bar, closed form expressions for the time dependent Green's function were obtained in [3] by inverting the Laplace 
transform directly. They are combinations of Heaviside functions similar to those in the the D'Alcmbert's solution 
for the entire line. However, the intermediate case, where the standing waves can form between the non-transparent 
boundary and the internal damper, was left untreated in [3], because neither modal nor direct Laplace inversion 
worked for it. We remedy the situation in Section [3] of this paper using our new method. 

Let us describe a physical picture underlying it. Intuitively, the D'Alembert's solution and its analogs can 
be interpreted as tracing the evolution of an initial <5-impulse as it splits into leftward and rightward traveling 
waves that undergo partial reflections at the dampers. The final formula represents a superposition of original and 
reflected traveling waves. This picture remains valid whether the standing waves can form or not, and provides a 
representation complementary to the modal one. It is available even in the critical cases, including the intermediate 
case mentioned above. Moreover, this superposition of traveling waves reduces to D'Alembert type solutions when 
they exist. Of course, the number of reflections may grow indefinitely as time goes on, but at any finite time only 
finitely many reflected waves contribute. Thus, one obtains a closed form solution as a sum with finite but time 
dependent number of terms. In fact, D'Alembert type solutions exist precisely when each wave can undergo only 
finitely many reflections in all of time. D'Alembert sums reduce to them in such cases, hence the name. 

It may seem from the above description that computations eventually become intractable as the number of 
reflections grows with time. This would be the case under the usual method of reflections [6l Sec 3.13.6], but luckily 
one can obtain all multiply reflected waves at once by a straightforward analytic manipulation. Moreover, for 
the parameter values close to critical ones the reflection coefficients are small, and waves that underwent multiple 
reflections are strongly attenuated. This means that taking into account only waves reflected a small number of 
times already gives a good approximation to the solution. 

The paper is organized as follows. In Section[5]we apply our method to the simple case of a vibrating bar with no 
internal damper. Even so, the closed form solution we obtain, Eq. (|16[) , appears to be new. In Section[3]we solve the 
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intermediate case left untreated in 0J, when the internal damper is present but exactly one of the bar's boundaries 
is transparent. In this case the denominator of the Green's function is still a sum of two exponents. Section 3] 
discusses some implementation issues for computing the vibratory response using the sums of traveling waves. In 
particular, some care is needed when taking the time derivatives since solution formulas involve discontinuous step 
functions. Computer simulations and comparisons to modal and FEM solutions follow the discussion. Finally, in 
Section [5] we introduce general D'Alembert sums when the denominator of the Green's function is a sum of many 
exponents. As an illustration, the method is applied to a vibrating bar again, but with no transparent boundaries. 



2 Bar with no internal damper 

In this section we show how to solve initial-boundary problem Eqs.([T]),([2]) by the method of D'Alembert sums in the 
simplest case, when no internal damper is present. Let £,[•] denote the Laplace transform over the time variable. 
Setting the external force and the initial data to zero and taking the Laplace transform we get the following boundary 
problem for U(x,s) :— C[u(x,t)]: 

s 2 s 

- T U(x,s) + 2h 3 -S(x-a)U(x,s)-U xx (x,s) = 0, (3) 



UJ0,s)-h 1 -U(0,s) = and UJL, s) + h 2 ~ U(L,s) = 0. (4) 
c c 

Since no internal damper is present h 3 = in this case. We will explicitly compute the Green's function G(x, £, s) 
for this problem and expand it into a series weighted by exponents with negative real powers to invert the Laplace 
transform. 

To this end, define ip(x, s), ?p(x,s) as solutions to ^U(x, s) — U xx (x,s) — satisfying only the left and the 
right boundary condition from Eq. Q respectively. This defines them up to a constant multiple and we make them 
unique by normalizing <yj(0, s) — 1 — tp(L, s). One easily finds that 

/ sx \ / sx \ 1 3I SK 

<p(x, s) = co S h(-j + hi sinh (— J = 2 [( X + h i) e ~ hi)e~—] (5) 
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^(a, s) = cosh^^— — ) + /i 2 sinh^ g ^ — ) = ^ 



a?) s(i — x) 

(1 + h2)e—*— + (1 - ft 2 )e ~ 



(6) 



By definition, the Green's function G(x, £, s) for system ([3j)— (|4j) satisfies 



-j G - G xx (x, s) = 5 (a: - 0, 



(7) 



G w (0,Z,s)-hi-G(0,t,s)=0 and G X (L, £, s) + h 2 - G(L, £, s) = 0. 
c c 



(8) 



As a function of x, G satisfies % U{x, s) — U xx (x, s) = for a; < £ and x > £. But, by construction, any solution to 
% Z7 — JX^ = satisfying the left (right) boundary condition must be a multiple of ip(x,s) (ip(x,s) ). Therefore, 
G(x,£, s) is equal to j4^(x, s) on [0, £) and Bip(x,s) on (£,£] with A and B independent of x. At x = £ it is 
continuous, but has a jump in the first derivative to produce 6(x — £) in Eq.(|7|). Namely, G x (£ , £, s) — G^^", £, s) = 
— 1 since G xx enters Eq.([7| with minus. Therefore, A,B can be found from the system 

A<p(Z,s)-B^,s) =0 
A V /(Z,8)-B1/{£,8) =1. 

Solving for them in the matrix form we get 



-WW] 




where W[y, 



is the Wronskian. For convenience, we denote 



A(s) := Wfv?,^] = (l + /ii/i 2 )sinh( — + (hi + h 2 ) coshf — 

s V c / V c 



(1 + /n)(l + ft 2 )e^ - (1 - hx)(l - h 2 )e 



(9) 
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so that explicitly A~c -} and B — c -^r7~T- We can write the Green's function in a simple form 
sA(s) sA(s) 



I ip(x,s)ip(C,s), X < £ 
GMs) = - W) l (10) 

\ip(x, s)<p(£,s), x > £_ . 

To solve the original problem we need to invert the Laplace transform and find T(x,£,t) :— £~ [G(x,£, s)]. In 
the modal approach G(x,£,s) is expanded into partial fractions ^ s _^ m , where p are the poles of G, then one has 
explicitly Cr 1 jrz^yt, = ( m -i)\ eP *' ^ or ^( x > £' ^ ms l ea ds to the expansion over the eigenmodes, standing waves, 
of the problem. Instead, we will expand G(x, £, s) into different functions, that still can be explicitly inverted. 
Our expansion is obtained using the special form of the denominator of G(x, £, s). Divide the numerator and the 
denominator of (JTUJl by (1 + hi)(l + h 2 )ei L . Then 

G(x,Z,a) = - — [ l _ hi)(1 _ h2) (11) 
1 (i+M(i+?'2) e c 

where we combined the x < £ and x > £ cases into 

inX ' 4 ' Sj " 2 S ^ + + + (l + h 2 ) C + (l + h 1 ){l + h 2 ) e J' 

When the real part Re(s) > is large enough we see from Eq. lfTT]) that the Green's function is a sum of a convergent 
geometric series 

( s) "^1(i + M(i + mJ 



e - N(x,£,s). (13) 



n=0 

Denote 0(x, £, t) := £7 1 [JV(a:, £, s)], then inverting (fT2|) we have 



' /,J ^Jf(ct-2L+(z + 0) + ^-^tf(ct-2L+|a;-£|) )■ ( 14 ) 



2i + /i 2 V i + ^i 

Recall that by the time shifting property of Laplace transforms, £~ 1 [e~ a,s i ;l ] = H(t — a)£ _1 [F](< — a), where F is 
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any function and 

I 1, t > 

(15) 




is the unit step Heaviside function. Applying this to Eq.(fT5|): 



n=0 



The n'th term of the sum is only non-zero for t > — n due to the presence of the Heaviside factors. In particular, for 
any fixed t only finitely many terms are non-zero, and their sum gives a closed form formula for the time dependent 
Green's function. Let [-J denote the floor function which returns the largest integer not exceeding its argument. 
Then we can replace the infinite sum by a finite D'Alembert sum 



[ct/2Lj 



n=0 



(i + /n)(i + /i 2 )y v ' v ' s ' c 



Similar expressions were obtained in [7] by a different method. When the attenuation factor fe^Hfe^j is small, 
only first few terms contribute significantly to the sum even for large t. Note that Ri := is exactly the reflection 
coefficient at the i'th boundary, the amplitude ratio of the reflected harmonic wave to the original one [5]. The 
attentuation factor in Eq(|16|) is R1R2, and it gives the amplitude decrease after a pair of reflections, one at each 
boundary. Note that it is < 1 by absolute value for non- negative hi. If one of the boundaries is transparent, say 
hi = 1, only the n = term in the sum survives and we get a D'Alembert type formula: 



r(x,£,t) = ! 



H(ct-\x-t\) + ±—!±H(ct-(x + t)) 
1 + hi 



(17) 



When also hi = 1 and both boundaries are transparent, it reduces to the D'Alembert's solution for the entire line. 
Of course, even Eq. (|16|) gives a closed form solution for any fixed time, but the number of terms required grows 
indefinitely as the time increases. 

Thus, we have a choice between expansions into standing waves (eigenmodes) and traveling waves. In the 
critical cases only the latter is possible, and it behaves nicely near the critical values as well. In contrast, the modal 
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expansion requires increasingly many terms for a good approximation near the critical values, and breaks down 
completely at those values. One can therefore say that the D'Alembert sum representation is complementary to 
the modal one from 0]. 



3 Bar with internal damper and transparent boundary 

It follows from the discussion in [4] that in the absense of internal damper either the Green's function has no poles 
at all, or the eigenmodes span the entire space. This dichotomy no longer holds in general. Intermediate cases arise 
when the internal damper is present, but one of the boundaries is transparent, and it is in these cases that the 
D'Alembert sums become indispensable. Standing waves can still form between the damper and the non-transparent 
boundary, but not enough of them exist to approximate arbitrary functions. Consequently, one can find neither a 
modal expansion nor a D'Alembert type formula for the solution. However, a D'Alembert sum representation exists 
in this case as well. 

As in Section [2J it will be convenient to express the boundary eigenvalue Green's function in terms of two 
particular solutions to Eq.Q. Denote by ip a (x,s) (ip a (x,s)) the solutions that satisfy the left (right) boundary 
condition only, and are normalized to be 1 at the corresponding boundary. Consider tp a first. Since the damper at 
x = a does not affect the equation on [0, a) we have f a (x, s) — tp(x, s) on this interval by definition of <p(x, s). For 

2 

x > a our (p a satisfies ^ U — U xx — again, but there must be a jump in its first derivative at a to accomodate the 
2/13^ S(x — a)U term in Eq.([3]). Along with continuity at a, we have </? a (a, s) — <ys(a, s) = and tp' a {a, s) — <p'(a, s) = 
2ha-ip(a, s). Since the difference ip a — <p also satisfies U — U xx = 0we get a Cauchy initial value problem for this 
function with initial point x — a. By inspection, (p a {x, s) — ip(x, s) — 2/13 tp(a, s) sinh (| (x — a)) for x > a. One can 
compute ip a analogously or notice that by symmetry it can be obtained from ip a by changing x to L — x, a to L — a, 
and hi to h 2 . With the help of the unit step Heaviside function H(-), Eq. (fT5l) . the x < a and x > a cases can be 
unified into 




(18) 




(19) 
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To further aid in notation, we define 



A a {s) =-~ W[<p a ,i> a ] = A(s) + 2h 3 (p(a, s)ip(a, s) , (20) 



where W[ip a , ip a ] is the Wronskian, and compute explicitly 



A (s) = (l + h 1 h2 + h 3 (h 1 + h 2 ))smh(— ) +(/ii+/i 2 + /i3(l+/ii/i2))coshf— 



V c 



+h 3 (h 2 - to) sinh( S(L fl) ) + h 3 (l - h x h 2 ) coshf S(L a) 



V C / V c 



(21) 



We will mostly use the exponential form of this expression 

A a (s) = i [(l + /i 1 )(l + / l2 )(l + ^ 3 )e* i -(l-/ii)(l-^)(l-/j 3 )e-* i 

+ (1 - to)(l + ha) to e^ L - 2Q ) + (1 + to)(l - h 2 ) h 3 e~^ L - 2a A . (22) 

Now we are ready to find the Green's function G{x, £, s). By definition, it satisfies 

2 

* G + 2h 3 - S(x -a)G- G xx (x, s) = 6(x - £), (23) 
c 

G x (0,e,s)-/ii^G(0,^s) =0 and G X (L, f, s) + G(£, £, a) = , (24) 

and has different analytic expressions depending on relative positions of x, a and £. Consider the case a < £ first. 
As a function of x, G satisfies Eq.Q for x < £ and x > £. Therefore, it is equal to Aip a (x, s) on [0, £) and Bijj(x, s) 
on (£, L] with A and £? independent of a;. At x — £ it is continuous, but has a jump in the first derivative to produce 
5{x - £) in Eq.(|2"5j). Namely, G x (£+,£,s) - G x (£ _ ,£,s) = -1 because G xx enters Eq.flU]) with minus. Therefore, 
A, B can be found from the system 

Asp a {^a)-B^%a) =0 
Va(f.*)--B^,«) =1 
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Solving for them we get in the matrix form 



-W[(pa,i>] 



-<Pa fa, 



But for £ > a Eq.fTSJ implies that s) = ^(£,s), so W[<p ,^ ] = ^[^a^] = -f A a (s) by Ea.(j2"0]l, Therefore, 



and B = c 



sA a (s) 



so that 



sA a (s) 



<p (x,s)^>(£,s), sc < £ 
V>(x,s)<p (£,s), a: > £. 



(25) 



Analogously, for a > £ we have 



sA Q (s) 



V>a(a;,sM£,s), sc > £ . 



(26) 



It will be convenient to rewrite G in a form that is both more explicit, and makes its symmetry G(x, £, s) = G(£, x, s) 
manifest. To this end, we introduce 



(p(x,s)ip(£,s), x < £ 
(p(^,s)ip(x,s), x > £ 



and compute 



9^{x,L s) = - [(1 + hi)(l + h 2 )ei( L -\*-m + (i _ + h 3 ) e 5(i-(«+0) 



+ (l + / ll )(l-^ 2 )e-t( L -^ + «» + (l-/i.i)(l-/i 2 )e-^ L - |a: -« l) . (27) 
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Analogously, let 



g a j,(x,Z,s) := 



5inh^f(x-a)j tpi^s), x<£ 
sinhff (£ - a)) V(z>*)> x > 6 



and g S ip(x,^, s) := 



sinhU (£ - a)) ^(af,s), x < £ 
sinh/f (i - a)) ¥>(£,s), x > £. 



Then 



<? s v(z, £, a) = J [(1 + Me 1 C £ -«-l*-«l) - (1 + h 2 )ei (L+a ^ x+ ^ 



+(1 - A a )e-5< £+0 -( a, + 6 ^ - (1 - h 2 )e-c( L - a -\ x -^] ■ (28) 



g sv (x, £,8) = - [(1 + h 1 )e^ x+ ^- a ^ - (1 + ^ 1 )ecC«-l=-«D +(l - /n)e-- c °-l*-*l5 - (1 - /*i)e-*« 1B+e 5-a)] _ ( 29 ) 



Since the g v ^, part is common to all arrangements of x, a and £ we have jointly 



G(x,£,s) 



sA a (s) 



9tpil:(x,€,s) + 2h 3 H(x - a)H(£ - a) (p(a,s) g s ^(x,(,s) 
+ 2h 3 H(a - x)H(a - f) ip(a, s) g sv (x, £, s) 



(30) 



Note that the last two terms are non-zero only when x and £ are on the same side of a. Therefore, whenever a 
separates x and £ the Green's function reduces to the first term. 

The denominator A a (s) is not a sum of two exponential terms as in Eq.®, and we can not apply the same 
simple trick in general (however, see Section [5]). But, if one of the boundaries is transparent, both the numerator 
and the denominator simplify significantly. If say h 2 — I, the above expressions reduce to 



A„(«) = (1 + M(l + ha) d L + (1 - h x ) h 3 e f( L - 2a ) = (1 + h x ){\ + h 3 ) < 



(1 — hi ) h 3 s „ 



(1 + fcl)(l + fc 3 ) 



<W(*>6 s ) = — [(! + hi)e-i^ + (1 - ht)e-^+0 



(31) 
(32) 
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2<p(a, s) g si ,(x, £, s) = [(1 + /u)e^ Q + (1 - h^e-^} 



e c 



e -t(l*-fl+o) _ e -t(*+f-a) 



(1 + /ii)e- = l x -«l - (1 + h^eri^+t- 2 ^ + (1 - ^ 1 ) e -td :r -«l +2a ) - (1 - h x )e--^ x 



(33) 



1 r 



(1 + hi)ei« x+ ®- a '> - (1 + /i 1 ) e iC«-l^-^D + (i _ /j^g-f («-!*-€!) _ (i _ /nje-f (O'+O-'O 



e<= 
~2" 



(1 + /ii)e-^( 2Q -( x+ «» - (1 + fci)e-5 1* -6 ! + (1 - /n)e--t Jta -l a '- € l5 - (1 - /u)e^ (a;+5) l . (34) 



Analogously to Eq.lfTTj) in Section [21 we can now represent the Green's function as 



G(x,£,a) = 



N(x,£,s) 



-, , (1-fri) h 3 -2^a ' 
1 (l+/ii)(l+/i 3 ) C 



(35) 



where (see Ea.([3"D|)) 



N(x,S,a) 



2(1 + h 3 )s 



V/>( x '£' s ) + h 3 H(x - a)H(£ - a) n sl p(x,£, s) + h 3 H(a - x)H(a - £) n sip (x, £, s) , (36) 



and 



^ipip{x j s) 

n s ^(x,£,s) 
n sip {x,£,s) 



2 • _ 
(l + /ii)et L 
2 ■ 2ip(a,s)g S Tp 

(l + /ii)et L 
2 ■ 2V>(q, s)ff sy 

(l + /ii)e* L 



1 - h x 
l + /ii 



-f(*+«) . 



,-f(z+£-2a) 



1 - /li 
l + /ll 



e _|(2a-(x+fl) _ e -t|x-e| , 

1 + 



(37) 



One easily verifies that whenever exponents in Eq. (|36|) have positive real parts they are multiplied by zero Heaviside 
factors. Thus, the inverse Laplace transform Q(x, £,t) :— Cj 1 [N(x, £, s)] is well-defined and given by 



2(l + /i 3 ) 



e^(x, e, t) + h 3 H(x - a)H(£ - a) 9^{x, £, t) + h 3 H(a - x)H{a - f) Sip (x, £, i) , (38) 
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where 



6^(x, £, t) = H(ct -\x-£\) + i-jA H(ct -(x + 0) ; 

1- // 



(s, £, t) = #(ci - |x - ^1) - JT(rf + 2a - (a + 0) + z r 1 [H{ct - 2a - \x - f |) - J3"(rf - (» + 0)] ; (39) 

1 + h\ 

6 stp (x, f , i) = H(ct -2a+{x + £)) - H(ct - \x - £|) + ^— Ji [ff(d - 2<z + |x - CI) - (ct - (x + £))] ■ 

1 + n.1 

Expanding G(a;, £, s) into a geometric series and inverting it termwise we obtain a D'Alembert sum as in Eq. (|16[) 
r(«, f , *) = ( TT-^v^t) ' Vet - 2na) 9 (*, t - *») , (40) 



n=0 



(1 + h t ){l + h 3 ) 



where [-J is the floor function returning the largest integer not exceeding its argument. The structure of T is 
identical to that of T from Eq. (|3l>)) with —h 3 /{l + h 3 ) playing the role of the reflection coefficient at the internal 
damper. By solving —h 3 /{\ + hs) = (1 — h 2 ) /(l + h 2 ) one finds that on [0, a] effect of the internal damper is similar 
to the effect of a boundary one placed at a with h 2 = 2/13 + 1, albeit not identical to it due to the difference in 

e(s,£,f). 



4 Implementation issues 

Computation of vibratory response via D'Alembert sums to presence of distributional terms in the Green's function. 
We address them in this section, and also compare D'Alembert sums to modal expansions and FEM solutions. The 
Laplace transform of the original system (P)-© is 

4 U(x, s) + 2h 3 -6(x - a) U(x, s) - U xx (x, s) = s + ^S+^Elll + 5(x ~ a) u(x, 0), (41) 

c c c c c 

U x (0,s) -hi-U(0,s) = u(0,0) and UJL, s) + h 2 - U(L, s) = — u(L, 0). (42) 

c c c c 
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In addition to integrals with the Green's function, the solution U(x, s) involves additional boundary terms due to 
the presence of spectral parameter in the boundary conditions, see [I]: 



U(x, s) = —G(x, 0, s) u(0, 0) + —G(x, L, s) u(L, 0) + 2—G(x, a, s) u(a, 0) 
c c c 



+ l\ G (x,t,s)^dt+ [ L GMa) *&°) + rt' a) dt. (43) 



Since T(x, £, t) := C s 1 [G(x, £, s)] properties of Laplace transform immediately imply 



c; 1 [s G(x, e, s)} = r t (x, s, t) + r(x, & o) 5{t) 

C^lGix, £, s) pg, a)] = [ T(x, U-r) p(& t) dr. 
Jo 



Inverting Eq. (|43p . the solution u{x,t) to the initial-boundary problem Eqs. ([!])- @ can be written in the form 



u(x, t) = - \h lU (0, 0) T(x, 0, t) + h 2 u(L, 0) T(x, L. t) + 2h 3 u(a, 0) T(x, a, t) 
c L 



(44) 



r t (x,z,t)u(£,o) + r(x,£,t)u(£,o) 



t r L 



T(x,Z,t-T)p{&,T)dtdfr 



Since r(x, £, t) is not a smooth function some care should be taken in computing its derivative T t , which is distri- 
butional. If we wish u(x,t) to represent the classical solution, one has to understand T t in Eq. (|4"4"]) as the classical 
part of this distributional derivative Sec3.12.2]. Note also that we dropped the distributional term T(x, £, 0) 6(t) 
when writing Eq. ([4~4"|) . It reflects the jump at t — 0, but is not a part of the classical solution u(x, t). Furthermore, 
consider the traveling wave expansion Eq. (|16[) for V . Differentiating the sum termwise and using the distributional 
product rule we get a sum of expressions like 



H t {ct-2nL)@(x,i, 



2nL 



H(ct-2nL)G t (x,£, 



2nL 
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The first term gives a sequence of delta functions. As with T(x, £, 0) 8(t), their contribution is purely distributional 
and is not a part of u(x,t). Therefore, for the purposes of computing the classical solution we have 

Equation (|38p for 0(x, £, t) itself involves Heaviside functions like H(ct — \x — £|), that need to be differentiated. 
Almost by definition H t (t) = S(t), where 8 is the Dirac's function, but one has to be careful with formal differenti- 
ation here. For instance, H(ct — \x — £|) = H(t — -\x — £|), but formally applying the chain rule on the right leads 
to an incorrect answer. Using distribution theory one can verify that [51 Sec 3.11.1]: 

H t (ct -\x- C|) = c8{ct -\x- C|) = c8{ct - x + £) + c:S(ct + x-£); 
H t (ct -2a±\x- C|) = c8(ct - 2a± \x-£\) = cH(±ct =p 2a) [8(ct - x + f) + <5(ct + x - £)] ■ 

Integration of the corresponding terms in Eq. (|44l) over £ then reduces to evaluating the initial displacement it(£) = 
0) at C = x — ct, x + ct, etc. For convenience of the reader, we give explicit formulas for convolutions of u(£) 
with time derivatives of functions from Eq.(|39|) : 



■ m(C) <K = cu{x — ctj + cu(x + ct) + c — u(ct — x) ; 



dt s v ' y ' " 1 + hx 

I H((,-a) dds ^,€,t) u (Q d £ = cH ( x _ ct _ a } u ( x _ ci ) + c H(x + ct- a)u{x + ct) - H(ct -x + a)u(ct -x + 2a) 



c - — ^-H(ct - 2a) H(x -ct + a)u(x -ct + 2a) + H(x + ct- 3a)u(x + ct - 2a) 

1 + All 



r°° QQ ( x £ t) 

/ H(a - C) — „ u(g) d£ = H(ct + x- a)u(2a -(x + ct))-cH(ct-x+ a)u(x -ct)-cH(a-(x + ct))u(x + ct) 
J~og ot 

+ c- — —H(2a-ct) \H(ct-x-a)u(x-ct+2a) + H(3a-(x+ct))u(x+ct-2a) . 
1 + hi L J 

As an illustration of this technique, we consider the response of the system to an initial Gaussian pulse centered 
at 0.25L with hi — 0.5 and h% — 0.7. We set L — 1.8 m, c = 1.5 m/s and a — 0.5L here and for all subsequent 
calculations. The right boundary is transparent, h% = 1, i.e. waves pass through it with no reflection, sec Fig(2] The 
intermediate damper affects a wave traveling to the right by reducing its hight in accordance with the parameter 
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h 3 . We wish to compare oaches. 




Figure 2: The response of the system for h x = 0.5, h 2 = 1, h 3 = 0.7 and a Gaussian impulse centered at 0.25L. 

Figure [3] depicts the response of the system at t — 1.5 s using the exact solution Eq.(Pf(J]) (left), and approximation 
errors for the modal approach |3j and FEM (right). The modal solution was calculated with h 2 = 0.999 since h% = 1 
is a critical value of the system and the modal approach breaks down at it [4]. This underscores the generality of 
D'Alembert sums as they work for critical values of the system as well. Furthermore, for t — 1.5 s one only needs 
to evaluate two terms in Eq. (j40"l) to obtain the exact response of the system due to the presence of the Heaviside 
factor. 
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Figure 3: The response of the system at t = 1.5 s with hi = 0.5 and /13 = 0.7 (a), and errors compared to modal 
approach (with h 2 = .999) and FEM (b). 



On the other hand, for the modal approach at least ten eigenfunctions are needed to get within 0.003 of the exact 
value of the response. For the FEM approach the system needs to be discretized with at least 180 elements for the 
value of the response to be within 0.0005 of the exact one. This illustrates the advantage of D'Alembert sums for 
small values of time. However, for larger times this advantage is lost since more and more terms would be needed 
in Eq. (|40[) . whereas modal approach and FEM require the same computational effort. 



5 General sums of traveling waves 

While constructing D'Alembert sums in the previous sections we used in an essential way the fact that the denom- 
inator of the Green's function is a sum of just two exponents. In this section we will show that the same approach 
applies when more than two exponents are present, albeit the sums become more cumbersome. Thus, the method 
can be applied to a large class of linear vibration problems. We also address analytic issues regarding termwise 
inversion of Laplace transforms in our context. 

Consider a Green's function of the form G(x, £, s) = ^ j — , where A(s) = J2k t =o a k e °' kS with real powers 

at- Let 5q be the largest power and set := cik/ao, := 5q ~ak- Dividing the numerator and the denominator 
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of G by aoe a ° s we have 



G(x,Z,s)= " { *>S' 3) _ aka . (46) 



This is exactly the form that we used in Eqs. (lll[) . (|35[) . only in both of them we had m = 1. Suppose further 
that \cfc 1 e- Soa F(x,Z,8)\ < C < oo for s with large real parts. For instance, if F itself is a linear combination of 
exponents with coefficients and powers depending on x and £, then this amounts to assuming that the real parts 
of their powers do not exceed 5o. This is always the case in initial-boundary problems regular in the sense of 
Tamarkin, see [8j. Now set q(s) :— Yl^—i a fc e " QfcS an( i own " niinfcafc, then for s with large enough real parts 
Re(s) one has \q(s)\ < max?; |a/ £ |e~ QminRc( ^ s ^ < 1. Therefore, the denominator of Eq. (|46|) can be expanded into a 
geometric series that converges absolutely and uniformly for Re(s) > u> with large enough u > 0: 

oo / m \ n 

n=0 \k=l ) 

oo I 

= E £ (-ir-T^ f af •••a^e-^+-+"-«" ^N(x,£,s), (47) 

^ ^ ni!---n m ! 

n— niH \-n m —n 

where we used the multinomial formula. Recall that the inverse Laplace transform amounts to integration along a 
vertical line in the complex plane also with large enough real part w : 



1 



r(x,£,t) = i^-. G(x^,s)e st ds. (48) 



Under our assumptions , \N(x,£, s)\ < £j and the summations in Eq. (|47p can be interchanged with the integral 
in Eq. (l48"]) for Re(s) > w. In other words, the series in Eq. (|47|) can be inverted termwise. Let 0(x, £, i) := 
Cj 1 [N(x,^,s)} as before, then 

m m 

C- 1 [ e -( n ^+-+ n ^ s N(x^,s)} = H(t -5>fc«fe) £™fe<*fc) ' 

fc=i fe=i 

Heaviside factors truncate the sum to finitely many terms for any given t > 0. Indeed, we have H\t— Y^k=i n k a k 
for t < JZic—i n k&k- Since Y)u—i n^a.^ > a m j n nk — a m in" the only non-zero terms correspond ton < t/a 
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Thus, a general D'Alembert sum generalizing Eqs. dTBl) . (p7f is a double sum 

[i/QminJ . m m 

T(x,U)= E E (-!)"— ;ar---a^H(t-Y,n k a k )e(x,U~Y, n ^)^ ( 49 ) 

n=0 niH hn m =n 1 m k=l k=l 

where the floor function |_-J returns the largest integer not exceeding its argument. Note that for Eq. (|4"9")) to make 
sense all a k must be real. This is an essential restriction on the type of problems admitting D'Alembert sum 
representation. Note also that the internal sum in Eg. (1491) contains potentially m n non-zero terms for every n, 
but for to = 1 it reduces to a single term. This is what makes D'Alembert sums far more attractive from the 
computational viewpoint when the Green's function has only two exponential terms in the denominator. Still, for 
small t evaluation of Eq. (|49[) remains practical even when m > 1. 

In the example of Section [3] before setting fi2 — 1 we had 5o = L/c, Si = (L — 2a) /c, S2 = — (L — 2a) /c and 
a3 = —L/c. Therefore, 777 = 3 and u\ — 2a/ c, 012 — 2{L — a)/c, = 2L/c and a m i n = min(a,L — a). Hence, for 
the exact answer at time t > one only needs to add up terms with n < 2m in(a L-a) ■ Moreover, from Eq.(f22j) 

= (1 -fei)fe 3 = (l-h 2 )h 3 = (l-feiXl-W-frO ^ 

ai (i + /7 1 )(i + /7 3 )' fl2 (i + h 2 )(i + h 3 y as (i + hl )(i + h 2 )(i + h a )- [ ' 

The function Q(x,^,t) has the same structure as in Eq. (|38|) and can be obtained in the same way. Expressions 
for 8^, 6 S ^ and 9 sip are more cumbersome in this case and we omit them here. We derived them explicitly for 
simulations below using a computer algebra system. 

To illustrate the theory, consider the response of the system with h\ = 0.5, h 2 — 0.7 and h 3 = 0.6 to an initial 
Gaussian pulse centered at 0.25L, Fig01 One observes that both boundaries are not transparent anymore (cf. 
FigEJ). We will compare the modal and FEM solutions to D'Alembert sums. Since the value of I12 is not critical 
here the modal approach can be applied exactly. Figure [5] depicts the response of the system using the exact solution 
Eq. (|4"9")) (a), and approximation errors for the modal approach [3]and FEM (b). As in Section[3J D'Alembert sums 
require the least computational effort. 

For larger times more terms need to be retained in Eqs. (flU)) . (l40l) and (|49]) to get an exact solution. However, 
when both boundaries are almost transparent one obtains an excellent approximation by retaining just first few of 
their terms due to strong attenuation of reflected waves. As a result, D'Alembert sums require less computational 
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effort than the modaf or FEM approaches even for large times. This can be best observed when the system is 
subjected to an external force to prevent vibrations from damping out quickly. 




Figure 4: The response of the system for h\ = 0.5, h 2 = 0.7 and /13 = 0.6, and a Gaussian impulse at 0.25L. 

Figure [6] shows the response of the system excited by an external harmonic force at 0.25L with parameters 
hi = hi = 0.9 and /13 = 0.6. The force per unit mass is of the form Q cos(uit)5(x — 0.25L), where Fq = 1 is 
a constant force per unit length and u) = 4 is the circular frequency. The left and right boundaries are almost 
transparent since the parameters hi and hi are close to 1. Consequently, very little reflection occurs at the 
boundaries, see Fig|51 
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Exact solution 2 terms at t= 1 .50 s 
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x ig~ 3 Exact solution 2 terms vs. Modal 10 eigenfunctions 




(a) 



0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.! 

x(m) 



(b) 




Figure 5: The response of the system for hi = 0.5, /12 = 0.7 and /13 = 0.6 at t = 1.5 s (a), and errors compared to 
modal approach and FEM (b). 



The response of the system at specific locations along the bar is tabulated in Tabled] The exact solution for time 
t = 10 s is u(x, 10) with N = [t/a m i n \ = 8 needed in the D'Alembert sum (|^]). while u e (x, 10) is the modal solution 
obtained using N e = 20 eigenfunctions. The last column shows approximation errors of the modal approach. One 
can see that to be accurate to three decimal places a large number of eigenfunctions is required. On the other hand, 
since the attenuation factor is small only a few first terms contribute significantly to the D'Alembert sum for all 
times. 
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x(m) \ °" 5 W 



Figure 6: The response of the system subjected to a harmonic force at 0.251/ with hi = hi = 0.9 and h§ = 0.6. 



X 


u(x, 10) N =8 


u e (x, 10)at c=20 


U(x, 10)at =8 - U e (x, 10) We =20 





.1092149996 


0.1096197171 


-0.4047175e-3 


.25 


.1051399774 


0.1048271564 


0.0003128210 


.5 


.07081830919 


0.06992067969 


0.00089762950 


.75 


0.6268305013e-l 


0.6295056781e-l 


-0.26751768e-3 


1 


0.03522116331 


0.3507059793e-l 


0.00015056538 


1.25 


0.3783383815e-2 


0.00387535258 


-0.91968765e-4 


1.5 


-0.2927453703e-l 


-0.2931169625e-l 


0.00003715922 


1.75 


-0.4979635515e-l 


-0.4983493334e-l 


0.00003857819 



Table 1: Exact and modal solutions compared, values of u are in meters. 



This is displayed in Tabled where approximations to the D'Alcmbcrt sum is truncated at first N — 0,1,2 terms 
in Eq. (|49p . Already with just one term in the sum (TV = 0), the displacements obtained are within two decimal 
places of the exact solution. Every additional term improves the accuracy by roughly one decimal place. Therefore, 
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for N = 1 we already obtain the accuracy of three decimal places that required 20 eigenfunctions under the modal 
approach. 



X 


u(x, 10)at =8 - u(x, 10)at =0 


u(x, W) N=S - u(x, 10)iv=i 


u(x, 10) N=S - u(x, 10)at =2 





-0.18118200c-2 


-0.1547189e-3 


0.0000040951 


.25 


0.0008789927 


-0.1706443e-3 


4.57e-8 


.5 


0.00248510275 


-0.13059818e-3 


-0.278129e-5 


.75 


0.00031824360 


-0.10003355e-3 


3.3294e-7 


1 


-0.146856055e-2 


-0.4258162e-4 


0.263416e-5 


1.25 


-0.1863464691e-2 


0.000009535038 


0.000002816701 


1.5 


-0.146038577c-2 


0.5756854e-4 


0.179306e-5 


1.75 


-0.43193244e-3 


0.00008094973 


1.56e-9 



Table 2: Difference between exact and truncated D'Alcmbcrt sums, values of u are in meters. 



6 Conclusions 

We described the method of D'Alembert sums for solving initial-boundary problems for partial differential equations 
with constant coefficients, and implemented it for vibrations of a bar with boundary and internal dampers. The 
method presents clear computational advantages over modal expansions and FEM for small evolution times and 
in critical cases, when the system of eigenmodes is incomplete. In fact, since the eigenmodes are not involved at 
all the method is not sensitive to phenomena related to them, for instance there is no difference in its application 
to self-adjoint or non self-adjoint problems. However, in practice D'Alcmbcrt sums are complementary to modal 
expansions being ineffective when the motion is dominated by few eigenmodes and most effective when they fail to 
form altogether. 

D'Alembert sums also hold theoretical interest providing a straightforward way to obtain closed form solutions to 
a number of vibratory problems. In particular, problems with multiple internal dampers and higher time derivatives 
in the boundary conditions can be considered. One essential limitation is that the characteristic determinant A(s) is 
a sum of real exponents, this is typically not the case for equations of order higher than two such as the biharmonic 
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equation for beams. 

Fully analytic treatment given here will not be possible for equations with variable coefficients, but the underlying 
idea of inverting the numerator and the denominator of the Green's function separately is more general. Whatever 
exact or approximate method is used to compute the Laplace transform it might be advantageous in some problems 
to expand the denominator into a D'Alembert type series for inversion. 
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